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Abstract: Preheating after inflation may lead to nonthermal phase transitions with 
symmetry restoration. These phase transitions may occur even if the total energy 
density of fluctuations produced during reheating is relatively small as compared 
with the vacuum energy in the state with restored symmetry As a result, in some 
inflationary models one encounters a secondary, nonthermal stage of inflation due to 
symmetry restoration after preheating. We review the theory of nonthermal phase 
transitions and make a prediction about the expansion factor during the secondary 
inflationary stage. We then present the results of lattice simulations which verify 
these predictions, and discuss possible implications of our results for the theory of 
formation of topological defects during nonthermal phase transitions. 
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1. Introduction 

The theory of cosmological phase transitions is usually associated with symmetry 
restoration due to high temperature effects and the subsequent symmetry breaking 
which occurs as the temperature decreases in an expanding universe |], §, ||. A 
particularly important version of this theory is the theory of first order cosmologi- 
cal phase transitions developed in ||. It served as a basis for the first versions of 
inflationary cosmology ||, as well as for the theory of electroweak baryogenesis ||. 

Recently it was pointed out that preheating after inflation J7j may rapidly pro- 
duce a large number of particles that for a long time remain in a state out of thermal 
equilibrium. These particles may lead to specific nonthermal cosmological phase 
transitions || [|. In some cases these phase transitions are first order [JT0|, |TT| ; they 



occur by the formation of bubbles of the phase with spontaneously broken symmetry 
inside the metastable symmetric phase. If the lifetime of the metastable state is 
large enough for the energy density of fluctuations to be diluted, one may encounter 
a short secondary stage of inflation after preheating ||. Such a secondary inflation 
stage, if it occurs late enough, could be important in solving the moduli and grav- 
itino problems. In this respect secondary "nonthermal" inflation due to preheating 
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may be an alternative to the "thermal inflation" JTJJ, suggested for solving these 
problems. 

In this paper we will briefly present the theory of such phase transitions and 
then give the results of numerical lattice simulations that directly demonstrated the 
possibility of such brief inflation. We will also discuss possible implications of our 
results for the theory of formation of topological defects during nonthermal phase 
transitions. A detailed description of numerical methods used in our work will be 
given in the appendix. 

2. Theory of the phase transition 

Consider a set of scalar fields with the potential 

V{^x) = \{^-v 2 f + 9 ^\\ (2.1) 

The inflaton field has a double-well potential and interacts with an iV-component 
scalar field x\ X 2 = Y,i=iXi- F° r simplicity, the field x is taken to be massless and 
without self-interaction. The fields couple minimally to gravity in a FRW universe 
with a scale factor a(t). 

The potential V(4>, x) has minima at = ±v, x — and a local maximum in 
the direction at = x = with curvature = —Xv 2 . The effective potential 
acquires corrections due to quantum and/or thermal fluctuations of the scalar fields 

III, 

= ^A(0 2 >0 2 + y <*V + y (4> 2 )X 2 + (2.2) 

where we have written only the leading terms depending on and x- The effective 
mass squared of the field is given by 

m 2 = -m 2 + 3A0 2 + 3A(0 2 > + g 2 ( X 2 }, (2.3) 

where m 2 = Xv 2 . Symmetry is restored, i.e. = becomes a stable equilibrium 
point, when the fluctuations (<f) 2 ), (x 2 ) become sufficiently large to make the effective 
mass squared positive at = 0. 

For example, one may consider matter in thermal equilibrium. Then, in the large 
temperature limit, one has ((f) 2 ) = (xj) = ff- The effective mass squared of the field 

mleff = -™ 2 + 3A0 2 + 3A(50 2 ) + g 2 ( X 2 } (2.4) 

is positive and symmetry is restored (i.e. = is the stable equilibrium point) for 
T > T c , where T 2 = 3 ^ 2 ^ g2 3> m 2 - At this temperature the energy density of the 
gas of ultrarelativistic particles is given by 

ir 2 _ 4 24m 4 7V(T c )7r 2 
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Here N(T) is the effective number of degrees of freedom at large temperature, which 
in realistic situations may vary from 10 2 to 10 3 . We will assume that Ng 2 ^> A, see 
below. For g A < 96J ^^ 7T A the thermal energy at the moment of the phase transition 
is greater than the vacuum energy density V(0) = which means that the phase 
transition does not involve a stage of inflation. 

In fact, the phase transition with symmetry breaking occurs not at T > T c , 
but somewhat earlier ||. To understand this effect let us compare the temperature 
T c ~ m/{yNg) and the mass m x = gcj) of the x particles in the minimum of the 
zero-temperature effective potential at = v — m /y/X. One can easily see that 
m x > T c for Ng A > A. This means that for Ng 4 3> A the temperature T c is 
insufficient to excite perturbations of the fields x% a ^ = v - As a result, these 
perturbations do not change the shape of the effective potential = v. Thus the 
potential at T slightly above T c has its old zero-temperature minimum at = v, as 
well as the temperature- induced minimum at = 0. Symmetry breaking occurs as 
a first-order phase transition due to formation of bubbles of the phase with « v 
at some temperature above T c when the minimum at = v becomes deeper than 
the minimum at = 0, and the probability of bubble formation becomes sufficiently 
large. A more detailed investigation in the case N = 1 shows that the phase transition 
is first order under a weaker condition g 3 ^> A ||. 

In the case Ng A > 10 2 A the phase transition occurs after a secondary stage 
of inflation. In this regime radiative corrections are important. They lead to the 
creation of a local minimum of V(0, x) at = even at zero temperature, and the 
phase transition occurs from a strongly supercooled state ||. That is why the first 
models of inflation required supercooling at the moment of the phase transition || . 

In supersymmetric theories one may have Ng 4 S> 10 2 A and still have a potential 
which is flat near the origin due to cancellation of quantum corrections of bosons and 
fermions [O. In such cases the thermal energy becomes smaller than the vacuum 
energy at T < To, where Tq = 2 ^. 2 m 2 f 2 . Then one may have a short stage of 
inflation which begins at T ~ To and ends at T = T c . During this time the universe 
may inflate by the factor 

* = £ ~ icr>(^) 1/4 . (2.6) 

Similar phase transitions may occur much more efficiently prior to thermaliza- 
tion, due to the anomalously large fluctuations (0 2 ) and (x 2 ) produced during pre- 
heating m U . These fluctuations can change the shape of the effective potential and 
lead to symmetry restoration. Afterwards, the universe expands, the values of (0 2 ) 
and (x 2 ) drop down, and the phase transition with symmetry breaking occurs. 

An interesting feature of nonthermal phase transitions is that they may occur 
even in theories where the usual thermal phase transitions do not happen. The 
main reason can be understood as follows. Suppose reheating occurs due to the 
decay of a scalar field with energy density p. If this energy is instantly thermalized, 
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then one obtains relativistic particles with energy density 0(T 4 ) which in the first 
approximation can be represented as p ~ E 2 ({(p 2 ) + (x 2 ))- Here E ~ T ~ p 1//4 is a 
typical energy of a particle in thermal equilibrium. After preheating, however, one 
has particles and x with much smaller energy but large occupation numbers. As 
a result, the same energy release may create much greater values of (0 2 ) and (x 2 ) 
than in the case of instant thermalization. This may lead to symmetry restoration 
after preheating even if the symmetry breaking occurs on the GUT scale, v ~ 10 16 
GeV § 1. 

The main conclusions of [§, |9| have been confirmed by detailed investigation 
using lattice simulations in [[II], [1^, |14, llf . One of the main results obtained in 



was that for sufficiently large g 2 nonthermal phase transitions are first order. 
They occur from a metastable vacuum at = due to the creation of bubbles with 
7^ 0. This result is very similar to the analogous result in the theory of thermal 



phase transitions ||. According to [10|, the necessary conditions for this transition 
to occur and to be of the first order can be formulated as follows: 

(i) At the time of the phase transition, the point = should be a local minimum 
of the effective potential. From (|2.3| ), we see that this means that Ng 2 (x 2 ) > At> 2 . 

(ii) At the same time, the typical momentum of x% particles should be smaller 
than gv . This is the condition of the existence of a potential barrier. Particles with 
momenta p < gv cannot penetrate the state with |0| « v, so they cannot change the 
shape of the effective potential at |0| ~ v. Therefore, if both conditions (i) and (ii) 
are satisfied, the effective potential has a local minimum at = and two degenerate 
minima at « ±v. 

(iii) Before the minima at ~ ±v become deeper than the minimum at = 0, the 
inflaton's zero mode should decay significantly, so that it performs small oscillations 
near = 0. Then, after the minimum at |0| ~ v becomes deeper than the minimum 
at = 0, fluctuations of drive the system over the potential barrier, creating an 
expanding bubble. 

The investigation performed in |L0| confirmed that for sufficiently large g 2 and 
N these conditions are indeed satisfied and the phase transition is first order. One 
may wonder whether for g 2 ^> A one may have a stage of inflation in the metastable 
vacuum = 0. 

Analytical estimates of Ref. [Q suggested that this is indeed the case, and the 
degree of this inflation for iV = 1 is expected to be 

-~(t) 1/4 - < 2 - 7 ' 

a v A ' 

/ 4 \ 1/4 

which is much greater than the number 10 ( y J in the thermal inflation scenario. 
One could also expect that the duration of inflation, just like the strength of the phase 
transition, increases if one considers N fields Xi with iV 3> 1. 
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However, the theory of preheating is extremely complicated, and there are some 
factors which could not be adequately taken into account in the simple estimates 
of ||. The most important factor is the effect of rescattering of particles produced 
during preheating This effect tends to shut down the resonant production of 
particles and thus shorten or prevent entirely the occurrence of a secondary stage of 
inflation. Thus the estimates above reflect the maximum degree of inflation possible 
for a given set of parameter values, but in practice the expansion factor will be 
somewhat smaller than these predictions. The only way to fully account for all the 
effects of backreaction and expansion is through numerical lattice simulations. In our 
paper we used a generalized version of the method of lattice simulations developed 
in [15, 10]. In the next section we will describe the basic features of our method and 
describe our main results. A detailed description of the lattice simulations will be 
given in the appendix. 

3. Simulation Results and Their Interpretation 

In our paper we will take A ~ 10~ 13 , which gives the proper magnitude of inflationary 
perturbations of density |p| [16|. We assume that g 2 ^> A, and consider v ~ 10 16 GeV, 



which corresponds to the GUT scale. A numerical investigation of preheating in the 
model ( |2.1| ) was first performed in |fT0| . The authors found a strongly first order 
phase transition. The strength of the phase transition increased with an increase of 
g 2 /X and of the number N of the fields x%- However, for the parameters of the model 
studied in ]^0| {g 2 /X ~ 200) there was no inflation during symmetry restoration. 
This is not unexpected because the estimates discussed above indicated that the 
expansion of the universe during the short stage of nonthermal inflation cannot be 
greater than (y) 1//4 . 

Keeping in mind that A in this model is extremely small, one would expect that 
in realistic versions of this model one may have g 2 / A as large as 10 10 , which could lead 
to a relatively long stage of inflation. However, for very large g 2 /X our analytical 
estimates are unreliable, and lattice simulations become extremely difficult: One 
needs to have enormously large lattices to keep both infrared and ultraviolet effects 
under control. 

To mimic the effects of large g 2 , we considered a large number of the fields Xi- 
We have performed simulations for g 2 /X = 800 and N = 19. With these parameters 
the strength of the phase transition became much greater, and there was a short 
stage of inflation prior to the phase transition. The details of our calculation and 
an explanation of our methods are given in the appendix. Here we only present the 
main results. 

The simulation showed that the oscillations of the inflaton field decreased until 
the field was trapped near zero. It remained there until the moment of the phase 
transition when it rapidly jumped to its symmetry breaking value, as shown in Fig. |I[ 
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The trapping of the field occurred because of the corrections to the effective 
potential induced by the particles and x produced during preheating, just like in 
the theory of high-temperature phase transitions. In our case, however, this effect 
has some unusual features. 

To first order in g 2 , the leading contribution to the equation of motion = —V 
is given by g 2 4>{x 2 ), where 

2 N Jn k k 2 dk 



Here ojf. = J k 2 + g 2 {4> 2 + (0 2 )) is the energy of Xi particles with momentum k and 
rifc is their occupation number; is the homogeneous component of the field. For 



I ,.?*'?.■. ■ (3-2) 



2 * I Jk 2 +g 2 m 



<C J (0 2 ), one has 



This quantity does not depend on 0; it can be evaluated using our lattice simulations 
when the field oscillates near = 0. It leads to the usual quadratic correction to the 
effective potential, see Eq. (|2.2[). This correction adequately describes the change of 
the shape of the effective potential for smaller than the amplitude of the oscillations 
of this field, because most of the time prior to the moment of the phase transition 
this amplitude is much smaller than J (0 2 ). 
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However, if we want to evaluate the effective potential at all values of |0| from 
to v, rather than for similar to the amplitude of the oscillations, then one should 
take into account that for sufficiently large |0| the term g\<f>\ becomes greater than 
9\J (4> 2 ) an d than the typical momentum k of particles Xi- 111 this case the main 
contribution to (x 2 ) is given by nonrelativistic particles with uo^ k, +g|0|, and one 
has 

<X 2 > « f 2 = ^ , (3.3) 

2vr 2 J g\(f)\ g\(j)\ 

where n x is the total density of all types of Xi particles. This implies that at large 
|0| the effective potential acquires a correction 

6V « g\cf>\n x . (3.4) 

Thus, instead of being quadratic or cubic in |0|, as one could expect from the analogy 



with the high-temperature theory || |ll| , the corrections to the effective potential at 
large \<j>\ are proportional to |0| 0. 

The combination of these two types of corrections to the effective potential 
(quadratic at small |0| and linear at large |0|) leads to the symmetry restoration 
that we have found in our lattice simulations. 

It is instructive to look in a more detailed way at the small region near the time 
of the phase transition. The first of the graphs in Fig. ||] shows the oscillations of 
the field soon before the phase transition, whereas the second one shows these 
oscillations soon afterwards. 




Figure 2: The spatial average of the inflation field as a function of time in the vicinity 
of the phase transition. The left figure shows the field just before the phase transition 
and at the moment of the transition. The field oscillates with an amplitude approaching 
lCrV The right fi gure shows the field after the phase transition, when it oscillates near 
the (time-dependent) position of the minimum of the effective potential at ~ v. Time is 
shown in Planck units. 
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First of all, one can see that just before the phase transition the field oscillates 
with an amplitude three orders of magnitude smaller than v, which is a clear sign 
of symmetry restoration. Another interesting feature is that the frequency of os- 
cillations does not vanish as we approach the phase transition, but remains nearly 
constant. Moreover, this frequency is only about two times smaller than the fre- 
quency of oscillations after the phase transition, which is equal to \/2m. Note that 
the frequency of the oscillations is determined by the effective mass of the scalar 
field, which is given by the curvature of the effective potential: m| = V" . This 
means that at the moment of the phase transition the effective potential has a deep 
minimum at = with curvature V" ~ +m 2 , i.e. the phase transition is strongly 
first order. Such phase transitions should occur due to the formation of bubbles 
containing nonvanishing field 0. 

Indeed, we have found that this transition occurred in a nearly spherical region 
of the lattice that quickly grew to encompass the entire space. The growth of this 
region of the new phase is shown in Fig. [| The nearly perfect sphericity of this 
region is an additional indication that the transition was strongly first-order. In 
comparison, the bubble observed in the lattice simulations of for g 2 / A ~ 200 was 
not exactly spherically symmetric. 

1282 . 18 




Figure 3: These plots show the region of space in which symmetry breaking has occurred 
at four successive times labeled by the value of the conformal time t pr . 

The first order phase transition and bubble formation seen in our simulations can 
be understood as a result of gradual accumulation of classical fluctuations 6cj)(t,x). 
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These fluctuations stochastically climb up from = towards the local maximum 
0* of the effective potential. Consider the regions in which 5<p(t,x) > 0*. If the 
probability of formation of such regions is small because they correspond to high 
peaks of the random field 0, then these regions will have a nearly spherical shape 
and can be represented by spherical surfaces of radius R* (bubbles). If the radius .R* 
is small, gradient terms will prevent the field inside the region from rolling down 
towards the global minimum at = v (subcritical bubble). If -R* is large enough, 
the gradient terms cannot push the field back to the metastable state = and 
the field inside the bubble rolls towards the global minimum, forming a bubble of 
ever increasing radius. This process can be described within the stochastic approach 



to tunneling proposed in [17]. Typically, the gradient terms cannot win over the 
potential energy terms if > OQm^ 1 !), where corresponds to the effective mass 
of the scalar field in the interior of the bubble. This provides an estimate for the 
initial size of the bubble R* ~ O^m^ 1 ]). 

The phase transition occurs from a state with energy density dominated by the 
vacuum energy density V(0). Figure |] shows the scale factor a as a function of time. 
The curvature becomes slightly positive at the time before the phase transition, 
which indicates a short stage of exponential growth of the universe. Because the 
curvature is hard to see in figure |] we have also plotted the second derivative a in 
figure |^. While the inflaton is trapped in the false vacuum state, a becomes positive, 
indicating a brief stage of inflation. 




2 -10 



3 -10" 



4 -10" 



5 -10" 



Figure 4: The scale factor a as a function of time. In the beginning a ~ \/t, which is a 
curve with negative curvature, but then at some stage it begins to turn upwards, indicating 
a short stage of inflation. 
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12 12 12 12 12 

u 1-10 2-10 3-10 4-10 5-10 

Figure 5: The second derivative of the scale factor, a. A universe dominated by ordinary 
matter (relativistic or nonrelativistic) will always have a < 0, whereas in an inflationary 
universe a > 0. We see that starting from the moment t ~ 2 x 10 12 (in Planck units) the 
universe experiences accelerated (inflationary) expansion. 



Another signature of inflation is an equation of state with negative pressure. Fig- 
ure H shows the parameter a = p/p, which becomes negative during the metastable 
phase. At the moment of the phase transition the universe becomes matter domi- 
nated and the pressure jumps to nearly 0. 



. 4 
0.2 



-0.2 
-0 . 4 
-0 . 6 
-0 . 8 
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D 12 12 12 12 12 

u 1-10 2-10 3-10 4-10 5-10 



Figure 6: The ratio of pressure to energy density p/p. (Values were time averaged over 
short time scales to make the plot smoother and more readable.) 
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From the beginning of this inflationary stage (roughly when the pressure becomes 
negative) to the moment of the phase transition the total expansion factor is 2.1. As 
expected this is of the same order but somewhat lower than the predicted maximum, 

/ 2 \ 1/4 

pa 5.3. We can thus conclude that it is possible to achieve inflation for 
parameters for which this would not have been possible in thermal equilibrium (g A <C 
A). In our simulation we have shown the occurrence of a very brief stage of inflation. 
This stage may be much longer for larger (realistic) values of g 2 /X. However, to 
check whether this is indeed the case one would need to perform a more detailed 
investigation on a lattice of a much greater size. 



4. Nonthermal phase transitions and production of topologi- 
cal defects 

In this section we would like to discuss possible implications of our investigation for 
the theory of production of topological defects after preheating [l(| [13|, [14 . 



The bubbles that appear after the phase transition can contain either positive 
or negative field, <fi = ±i>. If bubbles of either type are formed with comparable 
probability, then after the phase transition the universe becomes divided into nearly 
equal numbers of domains with = ±v, separated by domain walls. Such domain 
walls would lead to disastrous cosmological consequences, which would rule out the 
models where this may happen ||. || . 

In general, the number of bubbles with = +v may be much greater (or much 
smaller) than the number of bubbles with <fi = —v. Then the domain wall problem 
does not appear because the bubbles with = +v would rapidly eat all their com- 
petitors with <p = —v (or vice versa). This may happen, for example, if the moment 
of the bubble production is determined by the coherently oscillating scalar field <fi. 
In such a case, after oscillating a bit near the top of the effective potential, the field 
<p may wind up in the same minimum of the effective potential everywhere in the 
universe. 

To investigate the domain wall problem in our model one would need to repeat 
the calculation many times with slightly different initial conditions or to make them in 
a box of a much greater size that would allow one to see many bubbles simultaneously. 
Fortunately, the results obtained in our study may be sufficient to give an answer to 
this question without extremely large simulations. 

First of all, according to the stochastic approach [17]] to the theory of tunneling 



with bubble formation [T8|, [HJ, the bubbles of the field <p are created as a result of 
the accumulation of long-wavelength fluctuations of the scalar field with momenta k 
smaller than the typical mass scale associated with this field, see Section [3] In 
our case this mass scale is related to the frequency of oscillations of the scalar field 
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at the moment of the phase transition. At that moment the leading contribution to 
the fluctuations ((f) 2 ) is given by fluctuations with momenta much smaller than m^. 



We calculated the value of the long- wavelength component of y (0 2 ) and found that 
it is (approximately) of the same order as the amplitude of oscillations of the field 
at the moment of the phase transition. The existence of a first-order phase transi- 
tion suggests that the probability of bubble formation must have been exponentially 
suppressed during the metastable stage. Such suppression would only occur only if 
the amplitude of fluctuations required to form a bubble of the new phase was much 



larger than J (0 2 ) |T7]], which would in turn mean the required amplitude was much 
greater than the amplitude of oscillations of 0. This suggests that the probability of 
the bubble formation is almost entirely determined by the incoherent fluctuations of 
the field rather than by the small coherent oscillations of this field. Consequently, 
the probability of formation of bubbles containing = +v in the first approximation 
must be equal to the probability of formation of bubbles containing = —v. 

To make this statement more reliable one would need to estimate the amplitude 
of the long- wavelength fluctuations of the field in a more precise way, which would 
involve using lattices of a greater size. However, there is additional evidence suggest- 
ing that the number of bubbles with positive and negative must be approximately 
equal to each other. 

Indeed, as we have seen, the curvature of the effective potential remained ap- 
proximately constant during dozens of oscillations of the field prior to the moment 
of the phase transition. This suggests that the shape of the effective potential and, 
consequently, the probability of the tunneling, did not change much during a single 
oscillation. Therefore one may expect that, within a single oscillation, the probabil- 
ity of a bubble forming when the oscillating field was negative was approximately 
the same as the probability of the bubble forming when it was positive. 

If the number of bubbles with positive and negative is approximately equal to 
each other, the phase transition leads to the formation of dangerous domain walls, 
which rules out our model p0| . If correct, this is a rather important conclusion 



which shows that the investigation of nonthermal phase transition may rule out 
certain classes of inflationary models which otherwise would seem quite legitimate 

III- 

But this conclusion does not imply that all theories where the nonthermal phase 
transition is strongly first order are ruled out. For example, one may consider a model 
l[) with being not a real but a complex field, = ^(0i + i0 2 ), |0| 2 = §(0i + 0i)- 



Since the main contribution to the effective potential of the field in the theory Q2T 
is given not by the field(s) but by the fields \ii we expect that this generalization 
will not lead to a qualitative modification of our results. In particular, we expect 
that for sufficiently large N and g 2 / A the phase transition will be strongly first order 
and there will be a short stage of inflation after preheating. However, in the new 
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model we will have strings instead of domain walls. 

A similar model in the absence of interaction of the fields with the fields \ 



was studied in [13. 14]. It was argued that even in this case infinite strings may 
be formed. The theory of galaxy formation due to cosmic strings is currently out 
of favor, but it is certainly true that cosmic strings produced after inflation may 
add new interesting features to the standard theory of formation of the large-scale 
structure of the universe [pi] 

The possibility of strongly first order phase transitions induced by preheating in 
models with g 2 ^> A adds new evidence that infinite strings can be produced after 
nonthermal phase transitions. Indeed, infinite strings may not be produced if the 
direction in which the field falls from the point = at the moment of the phase 
transition is determined by the oscillations of the field 0. If, just as in the case 
discussed above, the amplitude of these oscillations are much smaller than <J (0 2 ) at 
the moment of the phase transition, then infinite strings are indeed formed. 



5. Conclusions 

The results of our lattice simulation confirm our expectations that preheating may 
lead to nonthermal phase transitions even in those theories where spontaneous sym- 
metry breaking occurs at the GUT scale, v ~ 10 16 GeV. Some time ago this question 
was intensely debated in the literature. Some authors claimed that nonthermal phase 
transitions induced by preheating are impossible, and the notion of the effective po- 
tential after preheating is useless. In our opinion, Figs. 1, 2 and 3 give a clear answer 
to this question. In particular, Fig. 1 shows that 90% of the time from the end of 
inflation to the moment of symmetry breaking the field oscillates about = with 
an amplitude much smaller than v. This could happen only because the corrections 
to the effective potential induced by particles and \ change the shape of V(4>) near 
= 0, turning its maximum into a deep local minimum. 

In some theories, this effect may lead to production of superheavy strings, which 
may have important cosmological implications for the theory of formation of the 
large scale structure of the universe. In some other theories, these phase transitions 
may lead to excessive production of monopoles and domain walls. This may rule out 
a broad class of otherwise acceptable inflationary models. 

In this paper we have shown that under certain conditions a nonthermal phase 
transition may lead to a short secondary stage of inflation. It would be interesting 
to study the possibility that a secondary stage of inflation induced by preheating 
could help solve the moduli and gravitino problems. The answer to this question will 
be strongly model-dependent because gravitinos can be produced by the oscillating 



scalar field even after the secondary inflation [22, E3|. Independently of all practical 



implications, the possibility of a secondary stage of inflation induced by preheating 
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seems very interesting because it clearly demonstrates the potential importance of 
nonperturbative effects in post-inflationary cosmology. 

The authors are grateful to Julian Borrill for very useful discussions. This work 
was supported by NSERC and CIAR and by NSF grant AST95-29-225. The work 
of G.F. and A.L. was also supported by NSF grant PHY-9870115, and the work of 
L.K. and A.L. by NATO Linkage Grant 975389. 



6. Appendix: The Numerical Calculations 

Our lattice program solves the classical equations of motion for a set of scalar fields 
in a Friedman-Robertson- Walker universe. These fields include an inflaton <p coupled 
to a set of matter fields Xi- The scale factor a is also solved for self-consistently. In 
this appendix we describe the exact equations being solved and the method used to 
solve them. 

Unless otherwise specified, all variables names refer to their bare physical values 
measured in Planck units. The quantities as they appear in the program have been 
rescaled in several ways, and these values will be indicated with pr, as in <p pr . The 
relations between the rescaled variables and their bare values are given below. Also, 
when we want to refer to a general property of the fields and x% we wm use a to 
indicate a generic scalar field. 

The simulations we discuss here used a grid of 128 3 points. The spacing of the 
points was chosen so as to ensure that the ultraviolet cutoff imposed by the grid 
was higher than all physically relevant momenta in the problem. To check this we 
monitored the power spectra |0fc| 2 and |xfc| 2 throughout the run. On a log plot 
you can see that there is a momentum above which the slope of the power spectrum 
decreases sharply (i.e. becomes more negative), and we set our grid spacing to ensure 
that this cutoff was well below the cutoff of the grid. A sample spectrum for x, taken 
from the late stages of preheating, is shown in figure [?]. The kink that can be seen 
at k « k max /2 persisted from preheating through the phase transition, and a similar 
feature can also be seen in the spectrum. 

The grid spacing we used was 

dx = 0.133^^- « 1.26 x 10* Mp 1 (6.1) 

measured in comoving coordinates with a = 1 at the beginning of the simulation. 
The total box size, L = 128dx, was equal to slightly more than 10 Hubble radii at 
the end of inflation. 

As a further check that the ultraviolet grid cutoff was not affecting the physics 
we did some smaller runs (i.e. with fewer fields) using these same parameters and 
another run using a grid of 64 3 points with twice as large a grid spacing. Both the 
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Figure 7: The spectrum |xfc| 2 as a function of momentum, shown in the late stages of 
preheating. The momenta k are shown in units of the Hubble constant at the end of 
inflation. 

64 3 and 128 3 runs had the same total box size and they showed essentially identical 
behavior for the fields, suggesting that the 128 3 grid has more than enough modes 
in the ultraviolet. 



6.1 The Field Equations 

We work in a FRW universe with metric g^ u = diag{l, —a 2 , —a 2 , —a 2 }. The equations 
of motion for the scalar fields and Xi derive from the potential 



V 



A 



(6.2) 



4 v J T 

where v , A, and g are constant parameters and x\ is understood to include a sum- 
mation over i. The results reported here used v = 7 x 10~ 4 , A = 9 x 10~ 14 , and 
g 2 = 800A, and took \ to be a 19-component field. This potential gives rise to the 
equations of motion 



a ■ 1 



+ 3-0--V^+(A(^-^)+<?V 







Xi + 3-Xi 



1 



-y'\ l + g 2 <p 2 Xl = v 



(6.3) 
(6.4) 



These equations can be simplified by the following variable redefinitions 

V = 1T& Xi'P r = v pr = -j^v; x pr = \/\(l> Q x] T pr = V\(f) J ^; ge = y (6.5) 



^0 



^0 
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where <fto is the initial (bare) value of the field <fi and all other parameters come from 
the original equations of motion. These then give the rescaled equations of motion 



'pr V pr pr +[(j) pr -a v pr + ge Xi, pr ~ — ) <l>w = (6.6) 



72 j. / j.2 J2..1 , „„ 2 

A.) 

a 



a" 



Xi,pr ~ ^ 2 prXi, P r + [ 9& <ftp T ~ ~ ) Xi,pr = (6.7) 



where derivatives are with respect to the rescaled time and distance variables defined 
above. Note that all first derivative terms have been eliminated, and the dependence 
on the coupling constants A and g is now only through the ratio g 2 / A, denoted here 
by ge. The value of A itself appears only in setting the initial conditions (described 
below). The rescalings involving O do not affect the equations of motion. 

6.2 The Scale Factor Equation 

Our simulations calculate a single scale factor at each time, neglecting metric fluc- 
tuations. The equation for the scale factor a is derived ^from the following two 
equations 

Att 

a— (p + 3p)a (6.8) 

3 

d\ 2 87T 

- ) =~o-P (6-9) 
a) 6 

where p and p are the total energy density and pressure, respectively. For simplicity 
we will first solve these for a single scalar field a. The energy density and pressure 
can be derived ^from the energy momentum tensor 

Tpv = cr,^,u - ^g^g al3 (^, a cr,i3 + g^V(a). (6.10) 

Assuming the fields are isotropic this equation can be solved and compared to the 
usual energy momentum tensor of matter 

T£ = diag{p } -p, -p, -p} (6.11) 

to give 



p=-a 2 + —\Va\ 2 + V(a) (6.12) 

p= l -a 2 -^\Va\ 2 -V{a). (6.13) 

Plugging these expressions into Eq. ( |6.8| ) and using Eq. ( |6.9| ) to eliminate the a 
term gives 

■d = -2- + 87ra(^-\Va\ 2 + V{a)) . (6.14) 
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Switching now to the fields (f) and \% an d to the potential in Eq. 

a = -4 + 8 ™ (3^ V< ^ 2 + l v ^l 2 ) + 1 <> 2 " 0" + ^x?) • (6-i5) 



Finally we rescale all variables according to Eq. (|6.5|) and take a spatial average 
(denoted by "<>") over the grid 



a 



a' 2 87T0Q 



(6.16) 



6.3 Initial Conditions 

Although the field equations are solved in configuration space with each lattice point 
representing a position in space, the initial conditions are set in momentum space and 
then Fourier transformed to give the initial values of the fields and their derivatives 
at each grid point. It is assumed that no significant particle production has occurred 
before the beginning of the program, so Minkowski space quantum fluctuations are 
used for setting the initial values of the modes. 1 

The calculation of these modes is given in detail in |15|]. The result is that each 
mode (Tfc has a probability distribution given by 

P(a k , a* k ) oc e- 2 "***** (6.17) 

where u k = \/k 2 + m 2 . Separating a k into a magnitude and a phase, the phase has 
a uniform probability distribution and the magnitude has a Rayleigh distribution 

P(\a k \) oc \a k \e- 2 ^\ (6.18) 

The momentum k is simply the Fourier transform variable and the mass is given by 

< = (6-19) 



Taking into account the finite size of the box, L, and the discretization of the 
spatial points with spacing dx and converting to the rescaled variables defined above 
yields for the initial magnitudes 

\<f>k, P r\ = Pr . 1/4 (6.20) 

sj2dx% r [k 2 pr + 3) 



lr This approximation is usually very good for fluctuations with momenta greater than H. For 
smaller momenta one can find a better approximation after performing the field quantization in 
curved space, see e.g. p3J, but in our case the corresponding corrections only slightly modify the 
final result. 
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times a Rayleigh distributed random number with standard deviation 1. Note that 
the program values of L, dx, and k are defined by the rescaling of x, and recall that 
ge = g 2 / A. Finally the zero mode, which appears as a value uniformly added to all 
grid points at the beginning of the calculation, is set to for \i an d 1 for (since 

4>pr = 0/0o)- 

In order to set the field derivatives it is necessary to know the time dependence of 
the vacuum fluctuations being considered. In Minkowski space this time dependence 
is given simply by the term e~ lujt which suggests 

cjfc = -iuj k a k . (6.22) 

Converting to rescaled values of time, mass, and momentum gives 

0fc,pr = i( t>k, P r\jklr + 3 ( 6 -23) 

Xik, P r = iXik, P r\/k 2 r + ge (6.24) 

The initial derivatives of the zero modes of and Xi are se t to 0. 

Note that having been eliminated from the equations of motion, the only place 
A shows up in the calculations at all aside from the ratio g 2 / A is in the y/\ term in 
the magnitude of the initial fluctuations. The coupling constant g appears nowhere 
except in g 2 /X. Meanwhile the initial value of the field 0, i.e. 0o, is used in the 
rescaling of the field and spacetime variables but it appears neither in the equations 
of motion for the fields nor in the amplitude of their initial fluctuations. In fact 
it appears in two places in the calculation. The first is in the evolution equation 
for the scale factor, Eq. (|6.16|) . The second way the initial value of enters the 
calculations is subtler. We said that the zero mode of <fi' pr is initially set to 0. At first 
this may seem like a poor choice since the field initially must be rolling towards 
0. In fact the beginning of the program is supposed to represent the end of inflation 
when the slow roll approximation is no longer valid, so why should 0' be set to 0? 
The answer comes from the use of the conformal field pr oc a0. As time goes on 
is decreasing but a is increasing, and there is a moment when these two balance 
and 4>' pr is momentarily 0. By setting (p' pr = the program automatically sets the 
beginning of the calculation at this moment. For the potential described here this 
moment occurs when « .35M P , so the initial conditions implicitly use the value 
of 0o- In the one place in the program where it does appear (in the scale factor 
evolution) this parameter is set to .35. This value is also useful in converting the 
output of the program to physically meaningful units. 
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6.4 The Calculational Method: Staggered Leapfrog 

The differential equations derived above are solved using a staggered leapfrog algo- 
rithm in which the field values and their time derivatives are calculated at alternating 
times separated by a half time-step. The time step is kept constant throughout the 
calculation. (Since the program uses conformal time the physical time elapsed dur- 
ing each time step changes as the program progresses.) This method is stable for 
second order equations involving no first time derivatives such as the field equations 
we use, see e.g. [p5 |. However, extra care must be taken in solving the equation 



for the scale factor, Eq. (|6.16| ), since this does contain a first derivative term a'. A 
naive calculation using the leapfrog algorithm as described above would mean that 
a" would be calculated at time r as a function of a(r) and a'{r — dr/2). There is a 
solution for this, although in practice the evolution of a is so slow and smooth that 
this problem makes no practical difference. We avoid this problem, though, by using 
the two following equations 2 

a\r + dr/2) ss a'(r - dr/2) + dr o"(r) (6.25) 

o'(r) « - (a'(r + dr/2) + a'(r - dr/2)) . (6.26) 
Solving these simultaneously, using Eq. (|6.16|) for a", gives 



a' {t + dr/2) « -a'(r - dr/2) + ^ ( -1 + a/1 + 2 fl/ ( r dr l 2 ) dr + L dr 2 \ ( 6 27 ) 

dr \ V o a J 

where 

/Mr)) = ^ (i (|VV| 2 + IVx^I 2 ) + \ - a 2 <) 2 + ipe ^ 



(6.28) 

Since Eq. ( 3.16 ) needs to be solved only once per time step this correction involves 



virtually no added computational time. 



2 G.F. would like to thank Julian Borrill, who suggested this solution. 
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